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Superconductors, ideally diamagnetic when in the Meissner state, can also exhibit paramagnetic 
behavior due to trapped magnetic flux. In the absence of pinning such paramagnetic response is 
weak, and ceases with increasing sample thickness. Here we show that in multiband superconductors 
paramagnetic response can be observed even in slab geometries, and can be far larger than any 
previous estimate - even multiply larger than the diamagnetic Meissner response for the same applied 
magnetic held. We link the appearance of this giant paramagnetic response to the broad crossover 
between conventional Type-I and Type-II superconductors, where Abrikosov vortices interact non- 
monotonically and multibody effects become important, causing unique flux conhgurations and their 
locking in the presence of surfaces. 


Introduction 


The diamagnetic Meissner effect is one of the hallmarks of superconductivity, where applied magnetic field is 
ideally screened out of the superconductor when cooled below the critical temperature T^. However, many field- 
cooled experiments on various materials over the past two decades detected a paramagnetic response, i.e. enhanced 
magnetic field inside the sample, usually referred to as paramagnetic Meissner effect (PME) or Wohlleben effect. The 
materials in question range from elementary ones such as Nb Qi, to much more complex high-Tc cuprates M. 
One proposed explanation for the enigmatic origin of PME in cuprates is based on the d-wave symmetry of the order 
parameter and the idea that tt junctions formed due to Josephson coupling between grain boundaries can result in 
spontaneous current loops with significant magnetic moments l3“ 3 ■ ^ much simpler and more general explanation 
is the compression and trapping of magnetic flux on cooling. Using this picture, Koshelev and Larkin [l^ calculated 
the magnitude of PME in thin stripes of conventional superconductors, and concluded that its theoretical maximum is 
~ 27% of the full Meissner response for the given magnetic field. Kostic et al. 17[ performed a supporting experiment 
on bulk Nb, and further concluded that polishing the sample surfaces strongly alters the PME, i.e. that surface 
barriers for flux entry and exit play an important role. 

The appearance of PME due to flux compression is easiest to understand in the case of mesoscopic samples, where 
the influence of confining boundaries is crucial. There, in analogy to surface superconductivity, during field-cooling 
the superconducting order parameter nucleates at the sample surface, and traps a multiquanta (giant) vortex inside 
the sample. Such large and compressed flux may lead to paramagnetic response, as first predicted by Moshchalkov et 
al. using self-consistent Ginzburg-Landau simulations |18l |. and subsequently verified experimentally by Geim et al. 
[H. A simple consideration shows that the paramagnetic moment strongly depends on the sample thickness, so that 
in very thick samples it appears only at very large fields and scales with penetration depth A over lateral size of the 
sample, while in thin plates it can be significant and scales with A over thickness 2^, 211. Therefore, the enigmatic 
PME in high-temperature superconductors becomes intrinsic for thin mesoscopic conventional superconductors. 

Recent years have seen the rise of interest in multiband superconductivity, particularly since its discovery in MgB 2 
and in iron-based materials 22h24|. The former has the highest of intermetallics; the latter are layered and 
high-temperature superconductors. To date, there have been no investigations of the paramagnetic response in these 
materials. Instead, a lot of attention has been paid to their rich vortex matter [^, and their possible classification 
outside the Type-I/Type-II dichotomy 26, 27| due to observed non-monotonic vortex interaction [2^. Early works on 
single-b and superconductors already discussed the broad crossover between conventional types of superconductivity 
where Abrikosov vortices exhibit long range attraction and penetration of vortices manifests as a large 


magnetization jump from the Meissner state to the mixed state (see e.g. Ref. [S^)- The lower bound of the crossover 
is given by the Hc{T) = Hg-iiT) line in the parametric space {He being the thermodynamic critical field and Hc 2 
the upper critical field) |32h37I| . below which textbook Type-I behavior takes place (for He > He 2 only Meissner 
state is thermodynamically stable, unless mesoscopic effects are strong, see Ref. [1^). The disappearance of the 
long-range vortex attraction marks the end of the crossover domain, and conventional Type-II behavior is recovered. 
This picture was recently detailed and extended to the multiband case in Ref. (^ . It is clear that non-monotonic 
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FIG. 1: Oblique view of the sample, the superconducting slab of width w, very long in other dimensions (indicated by dashed 
lines), in parallel magnetic field H. 


vortex interaction and other interplay effects between condensates in multiband superconductors are bound to also 
affect the interactions of trapped magnetic flux with the sample boundaries, and can lead to novel manifestations of 
the paramagnetic Meissner response. To reveal, quantify and explain the latter is the core objective of the present 
report. 


Results 


We consider a larger than mesoscopic two-band superconducting slab of width w {w/X ranges from 15 to 50 for 
the considered parameters) in a parallel magnetic field (see Fig. [T]), and report particular behavior of the sample 
magnetization as a function of the applied field [M(H) loops]. We primarily focus on two-band materials, but our 
findings can be qualitatively extrapolated to systems with more than two bands. The calculations are performed within 
the two-component Ginzburg-Landau (TCGL) theory (see Methods), where we have cautiously set a sufficiently high 
temperature T to ensure the qualitative and quantitative validity of our predictions in the context of recent debates 
and we used full microscopic expressions of all coefficients in the theory ^-4^. TGGL theory then comprises 
eight independent parameters, namely, the Fermi velocities of the bands vi and V 2 , the elements of the coupling matrix 
All, A 22 and A 12 = A 21 , the total density of states iV(0) as well as the partial density of states of the first band ni 
(note ni -I- ^2 = 1), and finally Tc, which sets the energy scale /7(^{3). By fixing the unit of length ^1 

and normalizing the order parameters by W, the parameters vi and Tc are fixed, and we are left with six parameters 
in the model: An, A 22 , A 12 , Vijv 2 , ni and fV(0). Instead of choosing fV(0), we opt to show the GL parameter of the 
first (stronger) band-condensate ki = '''^hich is an indicator for the expected magnetic behavior of the 

sample. 

In what follows, we consider an infinitely thick slab of width w = 8 OC 1 , and use periodic boundary conditions in 
the longitudinal direction (with size of the unit cell I = 120 ^ 1 , see Fig. [IJ. Without loss of generality, we take for the 
remaining microscopic parameters of the sample: ki = 1.5, An = 1.55, A 22 = 1.3, A 12 = 0.09 and m = 0.48. Note 
that such choice of parameters does not correspond to any particular material, and is actually by no means unique - 
since our main study will concern the dependence of the magnetic properties on the ratio of the Fermi velocities vi /f 2 
and temperature. We demonstrate these properties via calculated magnetization [M(H)] loops while adiabatically 
sweeping the magnetic held up and down. 

In Fig. [2]we show the M{H) loops at T = O.QdTc, for different values of vijv 2 - By increasing the latter parameter, 

we are actually decreasing the characteristic length scale of the second condensate C 2 = hv 2 /VQW (since Ci is hxed as 

2 

the unit of distance) and we are thereby driving the system into the Type-II magnetic behavior (since K 2 = -j- 


and the GL parameter of the coupled system k ^ are increasing; for calculation of the penetration depth A, please see 
Ref. [ 28 I) . This directly manifests in magnetization curves: for low vilv 2 {'^ 0.3) one easily recognizes typical response 
of a Type-I slab (see Fig. [DJa)), with superheated Meissner state in increasing held (with subsequent collapse to normal 
state), and supercooling in decreasing held [d^, with some hux trapping present; for high vi/v 2 (^ 0.65), the expected 
response of a Type-II slab is recovered [i^, still with some paramagnetic hux trapping (see Fig. [2](c)). However, at 
intermediate values of v\lv 2 a uniquely different shape of the magnetization loop is found, with a pronounced jump 
from the Meissner state to the mixed state with increasing held, and a very pronounced paramagnetic response in 
decreasing held (see Fig. [IKb)). 

In increasing held, all calculated magnetization loops exhibit a superheated Meissner state above the thermodynamic 





















3 



FIG. 2: Magnetization M{H) loops at T = 0.94Tc, for sequentially increased ratio of the Fermi velocities wi /«2 (and other 
parameters An = 1.55, A 22 = 1.3, A 12 = 0.09, ni = 0.48, and ki = 1.5), obtained by sweeping up and down the external 
magnetic field H (given in units of the thermodynamic critical field He). 


critical field He, where the superheating field Hsh agrees very well with the seminal calculations of Matricon and Saint 
James for Hsh{i^) of single-band materials [i^. At H = Hsh, superconductivity is either destroyed (for vi/v 2 < 0.34) 
or a jump to the mixed state occurs (for vi/v 2 > 0.34). The delimiting value of vi/v 2 = 0.34 exactly satishes the 
condition= Hc 2 - In decreasing magnetic field, the superconductivity nucleates at the surface superconductivity 
field Hc 3 |46|. Indeed, the nucleated states were only superconducting at the surfaces of the slab, with a large normal 
domain in the interior of the slab. For further lowered field and vijv 2 < 0.34 the normal domain remains trapped 
until abruptly expelled from the sample at the expulsion field H^- This analysis confirms that magnetic response 
of the system for vi/v 2 < 0.34 is the one of Type-I superconductors, since typical superheating-supercooling picture 
holds there, He 2 is smaller than He, and no vortices are found in the paramagnetic branch where flux was trapped 
upon nucleation of surface superconductivity. However, while decreasing held for ui/u 2 > 0.34, where consequently 
He 2 > He, the normal domain becomes unstable at held Hd but is not expelled; instead, it spreads into a vortex 
conhguration, stable down to persistently lower expulsion held He as vi/v 2 is increased. Simultaneously, hux trapping 
becomes notably more efficient, so that the vortex exit is hampered in decreasing held and paramagnetic response 
increases to its maximum at Hg. This tendency continues up to v\/v 2 ~ 0.53, for which paramagnetic response is 
almost an order of magnitude larger than the Meissner response a.t H = Hg, and approximately 30 times larger than 
the largest theoretical estimate of paramagnetic response to date (scaled to the diamagnetic response at a given held, 
see Ref. [l^). We therefore refer to this property as giant paramagnetic response (GPR). For vi/v 2 > 0.53, the 
cumulative paramagnetic response is still very large but gradually decreases, and magnetization curves in decreasing 
held connect to zero without any abrupt hux expulsion. In other words, we approach the Type-II limit, in which 


magnetization is expected to hover around zero for descending held in the presence of surface barriers |48| . In Fig. 
[21 we summarize the observed maximal amplitude, Max(47rM/i/) in the entire held range, and the total cumulative 
paramagnetic response, {AttM/H) = ^ Jh {M/H)dH, as a function of v\/v 2 , extracted from Fig. [21 

Based on Fig. |21 we argue that the giant paramagnetic response is characteristic for superconductors between 
conventional Type-I and Type-II (^ . Namely, this pronounced paramagnetic response is exactly found for sample 
parameters between the line HgiT) = He 2 {T) and the line where long-range vortex interaction changes sign (deter¬ 


mined by ehective GL parameter k* calculated after Ref. 28|), with a maximum found close to the parametric line 


where surface energy (asN) of the superconductor-normal metal (S-N) interface changes sign (determining the change 
in the polarity of the short-range vortex interaction (^). For the microscopic parameters considered here, we show 
this domain in Fig. |4l(a), as a function of vi/v 2 and temperature. To test our hypothesis further, we calculated an 
additional set of M{H) loops, shown in Fig. SKb), for fixed vi/v 2 = 0.55 and varied temperature indicated by yellow 
arrow in Fig. IlKa). From Fig. SKb), we confirmed exactly the same behavior of the loops and relationship of the giant 
paramagnetic response (GPR) with the delimiting lines of the critical domain: for T > O.QSTg the expected response 
of a Type-II slab is found, for 0.98Tc > T > O.dlTg the paramagnetic response in decreasing field increases when 
crossing the long-range vortex attraction line, and finally a pronounced paramagnetic response followed by a jump to 
the Meissner state is observed when crossing the 0 's at = 0 line. Besides being useful for reaffirming our conclusions. 
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FIG. 3: Maximal paramagnetic response in decreasing field at T — 0.94Tc (red) and its total cumulative value over the field 
span (black), as a function of Ui/ii 2 . Vertical lines indicate where Ha = Ha 2 , where the S-N surface energy changes sign (i.e. 
usN = 0), and where long-range interaction of vortices changes sign (left to right, respectively), delimiting the crossover range 
between standard types of superconductivity. 



FIG. 4: (a) The boundaries between different types of superconductivity in the {v\/v 2 ,T) plane, for other parameters as in 
Fig. m K* — l/%/2 line marks the onset of long-range attraction between vortices. At Ha{T) = Hc 2 {T) line, the mixed state 
vanishes in the bulk material. Dashed line shows where the energy of the superconductor-normal metal interface (osn) changes 
sign. Arrow shows the path to obtain the sequence of magnetization curves shown in panel (b), for ui/u 2 = 0.55 and varied 
temperature. Distinct changes in M{H) loops are found when either curve in panel (a) is crossed. 


this temperature dependence of the GPR can also be directly verifiable experimentally. Here, the considered samples 
are ideally clean, but even in realistic samples where flux trapping is present even at zero field, the rise and fall of GPR 
as a function of temperature will be easily observable in the above discussed scenario. Note that in general, changing 
any of the parameters can drive the in silico material across the crossover between the types of superconductivity, and 
thereby change the paramagnetic response. GPR is only sensitive on the regime of superconductivity the material is 
in, i.e., where the taken parameter set lies in the reconstructed Fig. |4l(a) - Type-I, Type-II superconductivity, or in 
between. 
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Discussion 

What is the underlying mechanism for the giant paramagnetic response? In simple terms, it is the facilitated 
trapping of magnetic flux in the crossover domain between Type-I and Type-II superconductivity, since vortices 
attract in the entire range of parameters where GPR is observed. However, GPR is found to be particularly large for 
ctsat > 0, where vortex-vortex interaction is purely attractive and vortices should coalesce into larger normal domains. 
On the contrary, we observe that in decreasing field separate vortex cores are still visible, though strongly overlapping 
(see inset in Fig. [5]). To clarify the dense vortex packing observed in Fig. [Sj we calculate the multibody vortex-vortex 
interaction shown for several vortex clusters in Fig. |6l As a major surprise, we found that in this regime multibody 
vortex interactions become short-range repulsive and cause the formation of a vortex lattice. This is illustrated in Fig. 
(HKa) (for vxjvi = 0.47 and T = O.OdTc, i-e. asN > 0), where we show the calculated vortex-vortex interaction as a 
function of the distance between vortices (labelled d), for a vortex pair, a vortex trimer, a vortex diamond-like cluster 
and a hexagonal vortex cluster. The pairwise vortex interaction is purely attractive, as expected, but in the other 
cases the short-range repulsion arises so that energetically favorable vortex-vortex distance arises in mid-range (note 
that this favorable distance closely corresponds to the average vortex distance observed in Fig. EKb)). An insight into 
the physics of this short-range repulsive interaction can be achieved by analysing the superconducting state inside, for 
example, the hexagonal vortex cluster shown in Fig [51 With this aim, we computed the maximum of the Cooper-pair 
density, Umax, inside that cluster for each band-condensate separately, shown as a function of vortex distance d in 
Fig. EKb). We reveal that the Cooper-pair density in the second condensate vanishes inside the vortex cluster at the 
vortex distance where short-range repulsion arises. Hence we can conclude that inside the vortex cluster the physics 
is driven by the other condensate, which has Type-H character, hence the repulsive interaction of vortices prevails 
at short distances. It is known that multibody vortex interactions are more complex than a simple superposition of 
pairwise interactions (see Refs. but it has never been found before that multibody interactions can change 

the polarity of the vortex-vortex interaction. This is a key feature of the found mixed state for parameters of the 
system between asN = 0 and Hc{T) = Hc 2 {T) lines in Fig. SI In addition, we have plotted in Fig. (Ua) the number 
of vortices in the sample Ny as a function of H in the downward branch of M{H) in Fig. [2] for vi/v 2 = 0.65 (in 
the Type-H limit) and vi/v 2 = 0.47 (inside the crossover region). The high retention of flux is clearly seen as a 
nonlinear behavior for v\/v 2 = 0.47 which contrasts the Type-H case in which Ny is linearly decreasing towards 
the origin. We find that although the number of vortices in the states between asj^ = 0 and Hc{T) = Hc 2 {T) lines 
slowly decreases with decreasing magnetic field, their favorable distance is approximately independent of field [see Fig. 
|5](b)]. This unconventional vortex state allows the penetration of the magnetic field in larger portions of the sample 
(inhomogeneous penetration, within but also between vortices), and clearly traps more flux than an ordinary vortex 
lattice, down to very low field - resulting in a more pronounced GPR. Due to the interlocking of vortices in this regime, 
the barrier for the expulsion of the entire vortex cluster in decreasing field corresponds to the Bean-Livingston barrier 
for a single vortex, which we confirmed by an independent calculation. Notice that as soon as the S-N surface energy 
changes sign, the barrier for single-vortex expulsion becomes nonzero at all fields. However, we have the simultaneous 
appearance of short-range vortex repulsion, which in effect diminishes the Bean-Livingston barrier and vortices are 
gradually expelled from the sample depending on their density and applied magnetic field. This manifests in the 
magnetization curves as a gradual decrease of the paramagnetic effect in decreasing field, down to zero for zero field. 
As the vi / V 2 ratio or temperature are further increased, vortices become increasingly repulsive and the paramagnetic 
response decreases to its conventional behavior for Type-H superconductors. 

In summary, we revealed a possibility of giant paramagnetic response in slabs of multiband superconductors (to 
which many recently discovered metal-borides, iron-chalcogenides, iron-pnictides, belong), with magnitude similar 
or multiply larger than the Meissner response for the same applied magnetic field. We showed that such unique 
magnetic response occurs in the crossover region between conventional types of superconductivity, and is not captured 
by the standard textbook descriptions. On technological end, our findings open a new class of desirable materials 
which can be switched to either strongly enhance or fully remove the applied magnetic field while having low power 
consumption. Further work is needed to characterize the behavior of these materials under e.g. applied electric current 
and nanostructuring or downscaling. 


Methods 

In this work we used the two-component Ginzburg-Landau (TGGL) theory. In the TCGL framework, as given in 
Ref. [ 2 ^ eight independent material parameters are needed for a system with both interband and magnetic coupling, 
namely, the Fermi velocity of the first band ui, the square of the ratio of the Fermi velocities in the two bands a = (^)^, 
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FIG. 5: The number of vortices in the sample in decreasing magnetic field (below H — Hd) at T = 0.94Tc (a), and the 
average distance between vortices (d^), for two values of ui/u 2 ratio that provide different sign of the supercondncting-normal 
state interface energy (asN)- Insets show cnmulative Cooper-pair density plots (|V'i|^ + IV' 2 p) of vortex states obtained in two 
considered cases for the same magnetic field H = 0.5G3Hc- 
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FIG. 6: (a) The vortex-vortex interaction energy, as a fnnction of the distance between vortices, for parameters leading to 
pairwise vortex attraction (crsjv > 0, see open dots). Nevertheless, the short-range repulsion between vortices arises for clusters 
comprising more than two vortices (insets depict the cnmulative Cooper-pair density distribution for the considered clnsters). 
(b) Maximum of the Cooper-pair density, Umax, inside the hexagonal vortex cluster for each band-condensate separately, shown 
as a function of vortex distance d between vortices. The shaded area delimits the short-range repulsion found for the hexagonal 
vortex cluster. 


the elements of the coupling matrix An, A 22 and A 12 = A 21 , the total density of states iV(0) as well as the partial 
density of states of the first band ni ( 71-2 = 1 — rii), and hnally T^, which sets the energy scale FF^ = 87r^T^/7C(3). 
The TCGL free energy functional reads 


i=i.2 ^ 



r(V'ri/’2 + V'ii/'2) + 


{h-Hf 

Stt 


( 1 ) 
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where j = 1,2 is the band index, aj = —N{0)njXj = —N{0)nj{T — Sj/rijS), Pj = {N{0)nj)/W"^, mj = 


3Wy{N{0)n,v]), and T = (iV( 0 )Ai 2 )/^, 
defined as in Ref. |4l|. The local magnetic 
Minimization of the free energy in Eq. 
for the order parameters ipi and and 
the coupled condensate). Introducing the 
^0 = hcfAe-KC,!, and for the lengths by Ci 


with 5 being the determinant of the coupling matrix, and S', Si and S 2 
field in the sample is denoted by h and the external applied field by H. 

(HD with respect to ipj and A yields the Ginzburg-Landau equations: two 
the equation for the vector potential (calculated from the supercurrent of 
normalization for the order parameters by W, for the vector potential by 
= hvi/\/6W, the dimensionless TCGL equations are written as: 


(-iV - Af^i - (xi - - 7V'2 = 0, 


1 




-(-iV -A) 4^2- (X 2 - \42\ )42 - 5 - 2^1 = 0, 




KjV X V X A = js 


where ki = 

hev^ 


yj 2 niN((A ^ ^2 = KiQ;-\/ni/n 2 , and 7 = Xi 2 /niS. In Eq. (|3]) the supercurrent density is 


I = - A)4i] + - A)r2], 


( 2 ) 

( 3 ) 

( 4 ) 

( 5 ) 


where TZ denotes the real part of the expression. After the made choice of normalization units, we are left with six 
parameters: An, A 22 , A 12 , Vijv 2 , Ui, and N(0). 

In our numerical experiment, the TGGL equations ©-(HI) were integrated self-consistently on a two dimensional grid 
with grid spacing ax = ay = much smaller than any characteristic length scale at the considered temperature. The 
discretization was implemented by the link variable method which preserves the gauge invariance of these equations 
M|. For the iterative solver, we combined a relaxation method with a stable and accurate semi-implicit algorithm 
55j. Periodic boundary conditions were applied in the x direction whereas for the y direction we imposed Neumann 
boundary conditions at the superconductor-vacuum interface (for details of the numerical implementation, please 
see Ref. [H^). Note that due to the infinite slab geometry, the surface magnetic field equals the applied one (the 
demagnetizing effects are negligible), and the simulation is effectively two-dimensional (in the (x, y) plane). The 
subsequently calculated magnetization, M = {{h) — H)/At: ((...) denotes spatial averaging inside the sample), is a 
measure of the expelled flux from the sample and the corresponding M(H) response was obtained by ramping up 
the magnetic field with steps of AH = 2 x 10“'^ (in units of Hq = hcj2e(^f). The magnetic field is scaled to the 
thermodynamic critical field for easier comprehension of the related physics. For calculation of He for two-band 
superconductors, we refer to Ref. [28 1. 
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